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ABSTRACT 

We analyze the dynamics of individual kilometer-size planetesimals in circumstellar orbits of 
a tight binary system. We include both the gravitational perturbations of the secondary star 
and a non-linear gas drag stemming from an eccentric gas disk with a finite precession rate. 
We consider several precession rates and eccentricities for the gas, and compare the results 
with a static disk in circular orbit. 

The disk precession introduces three main differences with respect to the classical static 
case: (i) The equilibrium secular solutions generated by the gas drag are no longer fixed points 
in the averaged system, but limit cycles with frequency equal to the precession rate of the gas. 
The amplitude of the cycle is inversely dependent on the body size, reaching negligible values 
for ^ 50 km size planetesimals. (ii) The maximum final eccentricity attainable by small bodies 
is restricted to the interval between the gas eccentricity and the forced eccentricity, and apsidal 
alignment is no longer guaranteed for planetesimals strongly coupled with the gas. (iii) The 
characteristic timescales of orbital decay and secular evolution decrease significantly with 
increasing precession rates, with values up to two orders of magnitude smaller than for static 
disks. 

Finally, we apply this analysis to the 7-Cephei system and estimate impact velocities for 
different size bodies and values of the gas eccentricity. For high disk eccentricities, we find 
that the disk precession decreases the velocity dispersion between different size planetesimals, 
thus contributing to accretional collisions in the outer parts of the disk. The opposite occurs 
for almost circular gas disks, where precession generates an increase in the relative velocities. 

Key words: planets and satellites: formation - stars: individual: 7-Cephei. 



1 INTRODUCTION 

The detection of giant planets in moderately close binary systems 
has raised many questions regarding the formation of these objects. 
For many years simulations of the dynamical evolution of circum- 
stellar disks suggested that planets may not form around the stars of 
a binary as the perturbation of the secondary star may (i) truncate 
the disk and remove the material that may be used in the forma- 
tion of planets, (ii) increase the relative velocities of planetesimals, 
which may cause their collisions to result in breakage and fragmen- 
tation, and (iii) destabilize the regions where the building blocks of 
these objects may exist. However, the discovery of planets around 
the primaries of the binaries 7 Cephei (Hatzes et al. 2003, Neuhuser 
et al. 2007), GL 86 (Els et al. 2001, Lagrange et al. 2006), HD 
41004 (Zucker et al. 2004; Raghavan et al. 2006), and HD 196885 
(Correia et al. 2008), where the stellar separation is smaller than 
20 AU, suggest that planet formation in such systems may be as 
efficient as around single stars. 

According to the core-accretion model (e.g. Safronov 1969, 
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Goldreich and Ward 1973, Pollack et al. 1996, Kokubo and Ida 
1998, Alibert et al. 2004, 2005), planetary cores are the results 
of accretional collisions between solid planetesimals immersed in 
the nebular gas disk. If this process is sufficiently fast to reach a 
few Earth-masses before the gas removal, gaseous envelops can 
collapse onto the embryos and result in giant planets. If growth 
is slower and/or if the bodies are located inside the ice line, the 
volatile component cannot participate in the accretion process, 
leading to the formation of small rocky terrestrial-type planets. 

Accretional collisions require low impact velocities in order to 
avoid disruption. In single star systems, relative velocities are usu- 
ally low, particularly during the early stages of planetary formation, 
and accretion appears as a natural outcome of most impacts. In bi- 
nary stellar systems, however, collisions are more complicated, es- 
pecially if the pericentric distance between stellar components is 
lower than ~ 20 AU. For instance, as the simulations show, it is 
possible to form terrestrial-class bodies from large ~ 1000 km- 
sized protoplanets around a star of close binaries (e.g. Quintana et 
al. 2002, Haghighipour and Raymond 2007). However, simulations 
of the collision and accretion of planetesimals have not been able to 
explain how these embryos formed in the first place (see Haghigh- 
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ipour 2008 for an extensive review). Tiie gravitational perturbations 
of the secondary star may affect the motion of bodies by excit- 
ing their orbits and increasing their eccentricities, which, for small 
kilometer-size planetesimals, may lead to encounter velocities well 
beyond the accretional limit (Heppenheimer 1978a, Whitmire et 
al. 1998). Had it not been for the fact that planets have already 
been detected in close binaries (e.g. 7 Cephei, GJ 86, HD41004A, 
HD196885), it would have been tempting to conclude that planet 
formation in these environments would be extremely unlikely. 

A possible mechanism for reducing the relative velocities 
of impacting planetesimals was presented by Marzari and Scholl 
(2000). After analyzing the effect of the nebular gas on the plan- 
etesimal dynamics, these authors found that the combined effect of 
the gravitational perturbation of the stellar companion and gas drag 
results in the appearance of an attractor in the dynamical system, 
causing the eccentricities of planetesimals to approach an equilib- 
rium value Beq. At this state, the pericenters of planetesimals orbits 
align and the difference between the longitude of the periastron of a 
planetesimal and that of the secondary star Accoq for all planetes- 
imals approaches equilibrium. Since encounter velocities depend 
on the dispersion in both the orbital eccentricity and longitude of 
periastron, at the equilibrium state the relative velocities of plan- 
etesimals are reduced and collisions occur below the fragmentation 
limit even (although the planetesimal swarm would still maintain an 
eccentric orbit). Subsequent studies by (Th6bault et al. 2004, 2008) 
show that the periastron alignment is size-dependent an is more ef- 
fective for planetesimals with equal sizes. The latter implies that 
even if at the onset of accretion, initial population of planetesimals 
consisted of equal-mass bodies, subsequent collision and growth of 
these objects would result in a mass spectrum which would lead to 
different values of ecq and AtUeq, causing the relative velocities of 
planetesimals to increase and bringing their accretion process to a 
halt. 

All the above-mentioned studies assume that the circumpri- 
mary gaseous disk is in circular motion. However, recent hydrody- 
namical simulations have portrayed a more complex picture. Apart 
from the truncation of the circumprimary disk near mean-motion 
resonances with the secondary star (ArtjTnowicz and Lubow 1994), 
and the subsequent alteration of the surface density, Kley et al. 
(2008) and Paardekooper et al. (2008) have shown that the gaseous 
disk develops an eccentricity Cg plus a retrograde precession with 
frequency Qg. The values of Cg and Qg appear very sensitive to the 
disk parameters (Kley et al. 2008). However, as long as the viscous 
time scale is much smaller than the induced precession time, the 
disk rotates as a rigid body and all fluid elements precess with the 
same rate. 

More recently, Marzari et al. (2009) found that self-gravity 
can damp the gas eccentricity and precession rate, and lead to an 
almost circular disk with no secular precession. However, these re- 
sults were obtained for disks of mass mo ~ 40Mjup and thus 
an order of magnitude more massive than those studied in Kley et 
al. (2008) or Paardekooper et al. (2008). Since Toomre's parameter 
is inversely proportional to mo, it is not yet clear if self-gravity 
would be relevant in disks with surface densities of the order of the 
MMSN. 

Paardekooper et al. (2008) studied the interactions of plan- 
etesimals in an static (i.e. no precession) eccentric disk in a binary 
stellar system. They found that, at any given semimajor axis, the 
value of the equilibrium eccentricity Coq depends on the size of the 
planetesimal. For small planetesimals (mainly coupled to the gas) 
this value is close to eg and for larger bodies it is close to the forced 
eccentricity e/. A similar behavior was also noted for Aweq. Inter- 



estingly, if Cg is sufficiently low, it is possible to find a critical semi- 
major axis for which = e/. In that case, all planetesimals would 
have the same equilibrium eccentricity independent of their sizes. 
More importantly, they would also exhibit apsidal alignment, con- 
stituting a very favorable breeding ground for planetary embryos. 

Unfortunately, two issues conspire against this idea. First, at 
least for the simulations done for the 7 Cephei system, the gaseous 
disk acquired too high eccentricity. Second and more importantly, 
the disk was assumed to be static. As can be seen from the hydrody- 
namical simulations by Kley and Nelson (2008), a precession in the 
disk appears to cause oscillations in the eccentricities of the plan- 
etesimals, which no longer seem to reach the stationary solutions. 

In this paper we aim to perform a detailed analysis of the dy- 
namics of individual planetesimals in a precessing eccentric gas 
disk that is also perturbed by the gravitational force of a secondary 
star. We will adopt the classical approach of considering a restricted 
planar three-body system in which the gas effects are mimicked 
with a non-linear drag force. Gas self-gravity will not be included, 
thus our results will only be applicable to light gas disks. Although 
such a model is a very simplified version of the real physical prob- 
lem, it has the advantage of isolating individual dynamical traits 
and allowing a separate analyze. As we will show, even this toy- 
model yields dynamical behavior quite distinct from that presented 
for circular gas disks. 

In Section 2, we construct the variational equations of the av- 
eraged secular system due to this drag. Although such expressions 
might have been used in an implicit form in hydrodynamical sim- 
ulations such as those mentioned above, we were unable to found 
explicit expressions for them in the literature. We present them in 
this paper, believing they would be beneficial to any reader wishing 
to imdertake similar studies. 

Section 3 is devoted to the evolution of individual orbits with- 
out an external perturber. We search for equilibrium solutions for 
a precessing disk and compare the results with the case of a static 
one. In Section 4, we include the gravitational perturbations of a 
binary stellar component and discuss its consequences on the dy- 
namics of the system. We show that the stationary orbits evolve 
towards limit cycles, in which both e and w oscillate aroimd equi- 
librium values with frequency equal to Qg. 

Finally, in Section 5, we apply this simplified analysis to the 7 
Cephei binary and estimate encounter velocities for different plan- 
etesimal pairs, as functions of semimajor axis and gas eccentricity. 
Conclusions and discussions end this articles in Section 6. 



2 GAS DRAG DUE TO AN ECCENTRIC 
CIRCUMPRIMARY DISK 

Our binary system consists of two stars with masses tua and free 
where niA is the primary. We place the origin of the coordinate 
system on ttia and assume that me has an orbit with a semimajor 
axis Ob and an eccentricity ee. We also assimie that both the gas 
disk and the disk of planetesimals orbit the primary star. Our focus 
will then be on the dynamics of planetesimals in the gaseous disk 
around the primary when subject to gas drag and the gravitational 
perturbation of the secondary star. 

The dynamics of the gas in its rotation arotmd the primary star 
is given by Euler's equation, 

dvg(r) ^ gruA ^ I dPgjr) 

dt Pg{^) dr 

In this equation, r represents the position vector of a gas element. 
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Pg(r) is the pressure of the gas, pg{r) is the gas density, Vg(r) is 
the velocity of the gas at position r, and Q is Newton's constant. 
The negative pressure gradient in equation lHJl causes each gas el- 
ement to orbit rriA with a sub-Keplerian velocity. Usually it is as- 
sumed that the rotation of gas is circular, where the deviation of the 
gas velocity from Keplerian (Vfe(r))circ is equal to 

~ - (2^)^ • 

The quantity r) = (AV(r)/Vfc(r))^;j.^ in such disks has been con- 
sidered to be between 0.005 and 0.01, suggesting that in a circu- 
larly rotating disk, the deviation from circular Keplerian veloc- 
ity due to pressure gradient is no more than 1%. In other words, 
Q = (1 + r;) ~ (Vg/V)£)circ > 0.99. Different authors have used 
different values for a. While Adachi et al. (1976), Gomes (1995), 
and Armitage (2010) considered ct — 0.995, Supulver and Lin 
(2000) mentioned that it is probably larger than ~ 0.99. 

In an eccentric disk, however, the situation is more compli- 
cated. A gas element in such disks rotates around the central star in 
an elliptical motion, and its Keplerian velocity varies with its po- 
sition vector. In other words, in eccentric gaseous disks we should 
adopt a — a(r). However, for the present simplified model we 
will restrict ourselves to a simple case with a constant value of 
Q = 0.995. 




X 



Figure 1. Astrocentric Cartesian coordinates showing the elliptic orbit of 
a solid planetesimal (black curve) immersed in an eccentric gas disk. The 
orbits of reference gas elements are shown by broad orange curves. 



2.1 Gas Drag 

The magnitude of gas drag is a function of the size of an object and 
its velocity relative to the gas Vrd ~ v — Vg. In this equation, v 
is the velocity of a planetesimal with a position vector r, and Vg 
is the velocity of gas at that location. For planetesimals in km-size 
range, gas drag is a non-linear function of the relative velocity and 
its magnitude is proportional to v^^i (Adachi et al 1976, Weiden- 
schilling 1977a, Supulver and Lin 2000, Haghighipour and Boss 
2003). The acceleration of a planetesimal due to the gas drag in 
this case is given by 



r = -C|vi.ci|v,.oi, 
where 

3Cd I V 



(3) 



(4) 



In equation lO, s is the radius of the planetesimal and pp is its 
volume density. The quantity Cd in this equation is the coefficient 
of the drag and has different functional forms for different values of 
the gas Reynolds number and planetesimal's radius. For km-sized 
planetesimals, Cd = 0.44 (Adachi et al 1976, Weidenschilling 
1977a). 



2.2 Relative Velocity 

We consider a gaseous disk with an eccentricity eg. We also as- 
sume that all gas elements have the same orbital orientation (i.e. 
longitude of pericenter vjg) at any given instant of time. Although 
hydrodynamical simulations (e.g. Marzari et al. 2009) have shown 
that both eg and vog can be a function of the radial distance, for our 
simplified model we have chosen constant values. This has been 
done primarily to avoid complications in the averaging of the vari- 
ational equations; however, it is fairly straightforward to implement 
a more realistic representation of the gas orbital dynamics. 

Figure[T]shows an schematic view of such a disk for eg = 0.2 
and an arbitrary value of the longitude of pericenter. The orbits of 



reference gas elements are shown by broad orange curves. The 
black ellipse represents the orbit of a solid planetesimal with a 
semimajor axis a and eccentricity e. The planetesimal's position 
is defined by its distance r from the primary star, its true anomaly 
/, and its longitude of the pericenter vj. 

To calculate the relative velocity of a planetesimal, we con- 
sider a two-body system consisting of the primary star and the plan- 
etesimal, and a Cartesian coordinate system with its origin at the 
location of the primary (Figure [T). In this system, the plane-polar 
coordinates of the velocity of the planetesimal are given by (Murray 
and Dermott 1999) 



Vr = r = 




(5) 



where yi = Qrap^ and p — a(l — e^) is the planetesimal's semi- 
lactus rectum. The mass of the planetesimal is assumed negligible 
with respect to ruA- 

As mentioned previously, we assume that the deviation of 
the velocity of a gas element from Keplerian is very small (a = 
0.995). As a first approximation, this assumption allows us to use 
equations similar to ^ to calculate the components of the velocity 
of the gas at the position of the planetesimal. As shown in Figure[T] 
the true anomaly of the orbit of a gas element at the radial distance 
r is related to / by 



fg ^ f + ZU — ZUg = f + AZU. 



(6) 



The plane-polar components of the velocity of gas at the position 
of the planetesimal can then be written as 



(7) 




1 + Cg COS (/ + Accj) 
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where Pg{r) — ag{r){l — e^). Given that in a two-body system 
r — a{l — e^) /(I + e cos /), the semimajor axis of the gas element 
can be written as 



1 - \ 1 + fig cos (/ + Azj) 



(9) 



, 1 + e cos / 

and subsequently 

1 + eg cos (/ + Aw) 
1 + e cos / 

The dependence on r is given implicitly through the true anomaly 
/. 

For the purpose of numerical simulations, it proves useful to 
calculate the contribution of the gas drag to the total acceleration of 
a planetesimal in Cartesian coordinate system. From equation {Sjl, 
we can write. 



= -C|Vrcl|(i- 

= -C|v,oii(j)- 



(10) 



To calculate the Cartesian components of Vg we note that from 
Figure [T] 



f + vj = fg + vjg = arctan(3//a;), 



(11) 



and therefore the Cartesian components of the gas velocity can be 
written as. 



— sin (/g + tUg) + eg sinci7g (12) 




cos (/g + ZUg) + eg COS ZUg 



In a precessing disk, the longitude of pericenter -cug is a function of 
time and can be included in numerical simulations as 



™3 = dgi + ^90' 



(13) 



where gg is the precessional frequency of the disk and zug^ is the 
value of longitude of pericenter at the beginning of the simulations. 



2.3 Averaged Variational Equations 

Finally, in this section we construct a simple analytical model for 
the differential equations averaged over the short period terms. In a 
perturbed two-body problem, the variational equations for (a, e, zu) 
can be obtained from Gauss' perturbation equations (Roy 2005): 

: [i?' esin/ + r' (1 + ecos/)] 



da 
di 
de 
It 



2a^ 



(14) 



^ \r' sin / + T' (cos / + cos it)] 

f - i/fhi?'cos/ + T'(l + ./p)sin/]. 

Here u is the eccentric anomaly of the planetesimal, R' is the radial 
component of the acceleration due to gas drag, and T' is its trans- 
verse component. From equation (O, these quantities are given by 

R' = -C\v^c;l\{Vr - Vrg) (15) 

T' = -C\Vrcl\{V0 - V0g). 

Substituting for the components of the planetesimal and gas veloc- 
ities from equations ^ and Q, expanding to the second order in e 
and Eg, retaining only lowest order terms in a and averaging over 
the planetesimal's mean anomaly, equations ( 114b can be written as 



for the semimajor axis, and 

dk ^ TT j ^ 

IE ^ ~ 
dh 



(16) 



= ~C^J^{k^kg)^/(k^kgfT(h'^h;^^ (17) 



dt 



C^j^{h-hg)^{k-kgY + {h-hgY., 



for the secular variables. Here we have adopted the regular vari- 
ables {k,h) — (e cos zu, e sin zu) for the motion of the planetesi- 
mal, and {kg,hg) — {eg coszug, eg sinzug) for the gas elements 
(Paardekooper et al. 2008). Due to its complicated form, the vari- 
ational equation for the semimajor axis J16b has been written only 
to its lowest order terms. We refer the reader to Gomes (1995) for 
more a detailed analysis of the complications of deducing an ana- 
lytical expression for valid for any value of a and eccentricity. 



3 DYNAMICS OF THE TWO-BODY PROBLEM 

Our first analysis is restricted to the dynamics of the planetesi- 
mal (around the primary star) without considering the perturbations 
from the binary companion. The orbital evolution of the planetes- 
imal will then be governed by the gravitational attraction of the 
primary star plus the effects of gas drag. In all the numerical sim- 
ulations performed in this section we assume a planetesimal with a 
bulk density of pp = 3 gr/cm'' and initial orbital elements of a = 2 
AU, e — 0.2 and Azu = 120°. The mass of the primary star is 
chosen as tha = 1.59M0, equal to the mass of the largest compo- 
nent of 7-Cephei. For the gas disk we assume a volumetric density 
of p = 5 X lO^^'' gr/cm'' (value at a = 2 AU), consistent with the 
values adopted by Paardekooper et al. (2008) for a MMSN. Finally, 
in our simulations both the gas eccentricity eg and the retrograde 
precession frequency gg are considered free parameters. 

Since we expect the motion of the planetesimal to be tightly 
coupled to the gas, we define a new set of regular variables 
{K, H) = {k — kg , h — hg) (Paardekooper et al. 2008) and re- 
write equations J17l ( as 

dK 



— = -c'kVk^ + H^ + gghg 
at 



(18) 



9gkg, 



dt 
where 

4 V a 



In the case of a static disk where gg — 0, equations ( I18t can 
be solved analytically to give 

KO rr,.^ Ho 



(19) 



Kit) 



1 + EoC t 



H{t)^ 



1 + EoC t 



(20) 



Here Ko and Ho are the values of K and H at t — 0, and Eo — 
Ko + Ho- The circularization time is thus proportional to 1/t, 
unlike the exponential decay time associated with the linear drag 
regime. The final equilibrium solution is given hy K = H = Q im- 
plying that, at equilibrium e — )■ eg and zj — >■ zog . Thus, the overall 
dynamical behavior of the planetesimal in a static eccentric disk is 
very similar to the circular case, except that the final equilibrium 
trajectory is now an ellipse. 

Neglecting terms proportional to the eccentricity, the rate of 
the change of the semimajor axis of the planetesimal can be written 
as 
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Figure 2. Top: Solution of the secular equations (18) for planetesimal with 
s = 1 km in a retrograde precessing disk with period of 1000 yrs. Middle: 
Orbital evolution, shown in the variables {k,h). Bottom: Same, but for new 
variables (K.H) = (e cos Aro, e sin Aro), where Azu = to — zug. In 
all cases the initial condition is marked with a filled black circle. 



= -C(l - a) , 



(21) 



dt 

where L — ^JJui is the Delaunay momentum associated to a. Since 
(1 — Of) <C 1, equation ( 121b suggests that the orbital decay of the 
planetesimal is much slower than the circularization time. In the 
limiting case where a = 1, this equations indicates that no sec- 
ular change exists in the semimajor axis of the planetesimal once 
the eccentricity and longitude of the pericenter have reached their 
equilibrium values. 



In the case of a precessing disk, system \\%\ is much more 
complicated to solve analytically. We can, however, focus on the 
equilibrium solutions. Figure |2] shows a typical example of the or- 
bital evolution of a one-kilometer planetesimal in a retrograde pre- 
cessing disk with period of 1000 yrs. The solution was obtained 
solving numerically the secular equations dlSt . Due to the disk 
precession, (K, H) no longer reach a stationary point but exhibit 
oscillations with frequency around a center located near the ori- 
gin. The same behavior is also noted for the time evolution of the 
original variables (fc, h) (middle plot). However, the system still 
retains an equilibrium solution in a new set of variables defined as 
{K,H) — (e cos Atu, e sin A-ro), where Aro = — zug. This 
behavior has been noted for all values of s and as well as for any 
initial condition. 

To obtain analytical expressions for the equilibrium solutions 
in this case, we note that from the definitions of (k, h) and (kg, hg), 
it is possible to write 

(22) 



The corresponding variational equations of K and H are then given 
by 



dK 
dH 



dk 
'dt 
dh 

'dt''' 



dh 
'dt 
dk, 
'dt'"' 



ggCgH 

QgCgK . 



(23) 



In an equilibrium state, the values of {K, H) are constant and equa- 
tions (HI can be simplify to 



dk 
dt 
dh 

'dt''' 



dh, 
'''+m'' 



dk, 

l-t''' ^ 



ggCgH 
ggCgK . 



(24) 



Equations i24\ are a set of algebraic equations that can be solved 
analytically. Substituting for the time derivatives of k and h from 
equations J17t . the equilibrium values of K and H are given by 



(25) 



Keg = 


e2 


eg 


Heq = 


^eq ^ eq 


where 





2,1 

"""^ 2\C 



9a 



9g 



4e2 



9s 
2C' 



(26) 



For a static gas disk (i.e. gg — 0), this equation indicates that 
the planetesimal's eccentricity e^q — > eg, whereas for an eccentric 
disk Esq varies according to the ratio between the disk precession 
frequency and the drag coefficient (i.e. object size). Figure[3]shows 
the equilibrium values of the eccentricity and Azu for different val- 
ues of planetesimal radius and three values of the gg. Even small 
values of the precession frequencies cause significant changes in 
the dynamics of the system. Except for very small sizes where the 
coupling to the gas is strong, the planetesimal does not follow the 
gas exactly. At large sizes, even though the eccentricity still reaches 
an equilibrium, its value is smaller than eg, reaching e^q — >■ for 
large bodies or faster precession rates. Similarly, even though the 
longitude of pericenter of the planetesimal is still locked to the pre- 
cession of the disk and circulates with the same frequency, the ap- 
sides are no longer aligned. After a transition time, Azu acquires 
an equilibrium value of Azu > (for retrograde precessing disks), 



6 C. Beauge, A.M. Leiva, N. Haghighipour and J. Correa Otto 







g =0 






\ \ 27i/lg^l=4000 yrs 






\l()()() \ 






\200 \ \ 











0.1 



1 10 
planetesimal radius [km] 



100 




1 10 
planetesimal radius [km] 



100 



Figure 3. Top: Equilibrium eccentricity, as function of the planetesimal 
radius s, in a disk with eccentricity eg =0.1 with different values of the 
precession frequency gg. Bottom: Equilibrium value of Aro = ■cu — TUg. 
We assume a retrograde precession (i.e. gg < 0). 



truncated up to second-order in the eccentricity. The complete av- 
eraged equations of motion for (k, h) are given by 

dk 



dk 



(28) 



dh 



dh. 



where the first terms are the contributions of the gas drag (expres- 
sions |[17)) and 



3 rriB niA + mB 

4 mA V 

5 acB 



ruA 



1)3/2 



4aB(l- 



(29) 
(30) 



In the absence of gas drag, the semimajor axis is constant and 
the solution of (fc, h) are given by the classical linear Lagrange- 
Laplace model 

p,i(gt + <#.o) 



k(t) + ih{t) 



+ e/, 



(31) 



where we have used the notation = expa;. In this equation 
and ef are the proper (or free) and forced eccentricities, re- 
spectively, and <^o is the initial phase angle. The expression for 
the forced eccentricity is given by ( I30t . and the free eccentricity 
is equal to 



ep = (ko - e/)2 + 



"•0- 



(32) 



Here ko, ho denote the initial values of their corresponding quan- 
tities at t = 0. Note that g is the secular frequency of the system. 
We refer the reader to Heppenheimer (1978b), Marzari and SchoU 
(2000), Thebault et al. (2004, 2006), and Paardekooper et al. (2008) 
for more details. For any initial values of the quantities k, h, and 0, 
the minimum and maximum values of the eccentricity are equal to 
emin = 1 6/ — Gp | and Cmax = 6/ + fip, respectively. In the case of 
initially circular orbits, Cmin = and emax = 2e/. 



indicating that the planetesimal always is behind the motion of the 
gafl 



4 SECULAR DYNAMICS IN THE RESTRICTED 
THREE-BODY PROBLEM 

4.1 Equations of Motion 

We now analyze the dynamics of the planetesimal including the 
gravitational perturbations from the stellar companion m b ■ Similar 
to the previous case, both the gas disk and planetesimal orbit the 
primary star tua, and the coordinate system will also be centered 
on this body. 

Outside mean-motion resonances, the dynamics is dominated 
by secular perturbations. Averaged over short-period variables (i.e., 
the mean anomalies), the disturbing function can be written in 
terms of regular variables (Heppenheimer 1978a), 



GmB 



8(1 



)3/2 



2aB(l 



(27) 



^ Direct precession of the disk implies Aro < (planetesimal trails disk) 
while retrograde precession implies Aro > (planetesimal precedes disk). 



4.2 Combined Effects with a Precessing Disk 

The complete dynamical system contains perturbations with two 
distinct frequencies: Qg and g. Except for large semimajor axes, 
we expect g <C \gg\ and resonances between the two frequen- 
cies should not to be significant. In order to have an insight on 
this complex problem, we show to examples of numerical integra- 
tions in Figure|4] Left frames correspond to a planetesimal of radius 
s — 0.5 km, while the right plots are drawn for s = 50 km. The 
figure caption gives initial conditions for the planetesimals and the 
parameters of the disk. Red symbols correspond to the solutions of 
the secular system ( |28t , while results from full N-body simulations 
are shown in black. The agreement between these two results show 
that the analytical equations give a very good representation of the 
real dynamics of the system. 

At first sight, the orbital evolution of both planetesimals ap- 
pear significantly different. For the smaller body, neither e nor Auj 
reach zero-amplitude equilibrium values, but evolve towards a limit 
cycle displaying significant oscillations with frequency gg. Never- 
theless, Azu librates while the longitude of pericenter of the plan- 
etesimal circulates. For the larger object the plots seem to indicate 
an opposite behavior. In this case, zu oscillates with a small ampli- 
tude around an equilibrium value, even if this behavior implies a 
circulation with respect to the gas disk. 

Figure [5] shows the solutions of equations l |28l l in the {k,h) 
plane, after the transient period is over and the orbit reaches the 
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Figure 4. Simulation of the orbital evolution of two planetesimals (s = 0.5 
km on the left, and s = 50 km on the right) under the combined effects 
of gas drag from an eccentric disk (eg = 0.1), precession of the disk 
(27r/|gg I = 1000 yrs) and the gravitational perturbations from a secondary 
star equal to that of 7-Cephei. Initial conditions for the planetesimals were 
a = 2 AU, e = 0.1 and Aro = 120°. Red lines show results from the 
analytical differential equations j28t while black corresponds to exact nu- 
merical integrations. Gas volumetric density at 2 AU was chosen equal to 
p = 5 X 10"^" gr/cm-^. 



limit cycle. Graphs liave been made for different values of the plan- 
etesimals radius s. There appears to be a smooth transition between 
small and large planetesimals with no topological changes in the so- 
lutions. The only difference is in the relative magnitude between the 
forced and free eccentricities. For small planetesimals, the limit cy- 
cle includes the origin of the coordinates and uj circulates, whereas 
for large bodies the opposite occurs and uj librates. 



5 RELATIVE VELOCITIES BETWEEN 
PLANETESIMALS 

5.1 Orbits of Different-Size Bodies 

In order to determine the effect of the precession of the gaseous 
disk on the process of the accretion of planetesimals, we will ana- 
lyze the relative motions of different size bodies in crossing orbits. 
Figure [6] shows the evolution of the orbital eccentricity for plan- 
etesimals with s between 0.5 and 50 kilometers. In each case we 
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Figure 5. Limit cycles in the {k, h) plane for various planetesimal radii. 
Initial conditions and gas parameters are the same as in the Figure|4] 
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Figure 6. Equilibrium solutions for static (left) and precessing disks (right) 
for two values of the gas eccentiicity eg . Semimajor axis of planetesimals 
was a = 2 AU giving a forced eccentricity of e ^ = 0.057. The gas density 
at that point was chosen to be p = 5 X 10^^" gxicw? . 



have only plotted the orbit after the limit cycle has been achieved. 
The top panels show results for a quasi-circular gaseous disk with 
eg = 0.02, while in the bottom graphs we chose eg = 0.2. For 
comparison, on the left-hand side we have plotted the results for a 
static disk (gg = 0) while on the right we chose a retrograde disk 
with a precessional period of 1000 years. In all cases the initial 
value of the semimajor axis of the planetesimal was set to a = 2 
AU which gives a forced eccentricity equal to e/ = 0.057. 

In a static disk (upper left panel), the equilibrium eccentrici- 
ties lie in the interval between eg and e / depending on the size of 
the planetesimal (Paardekooper et al. 2008). Since in the Simula- 
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tions depicted by the top panels the value of Cg is close to e/, there 
is little dispersion in the eccentricities and all orbits have similar 
trajectories independent of their sizes. However, for the more ec- 
centric disk (bottom panels), the spread in eccentricities is much 
more pronounced. Such a separation between objects of different 
sizes results in a larger relative velocities during physical encoun- 
ters. 

The effect of the disk precession can be seen comparing the 
right-hand plots. A non-zero value of gg lowers the equilibrium ec- 
centricity of a planetesimal to values below eg, even when the size 
of the planetesimal is small (see also Figure[3]l. In the case where 
the forced eccentricity is larger than eg, the precession of the disk 
increases the range of the equilibrium eccentricity of a planetesimal 
compared to its corresponding value in the static disk. This implies 
that for Sg 7^ the relative velocities of planetesimals would be 
larger than in a static disk. Conversely, if eg > e/ , such as depicted 
in the bottom panels, the precession will induce a smaller range in 
eccentricity for different size bodies. This can be seen by compar- 
ing the eccentricities of two planetesimals with sizes of 5 km and 
50 km in both the static and precessing disks. 

A second effect of the disk precession, when coupled with the 
gravitational perturbation from the binary companion (ms), is the 
appearance of limit cycles, causing the periodic time variations of 
planetesimals' final eccentricities. Fortunately, however, the phase 
angles of the limit cycles are only weakly dependent on the radius 
of a planetesimal s, which causes these oscillations to appear al- 
most coherent. 

In short, in disks with high values of eg, the disk precession 
can actually reduce the relative velocities of planetesimals with dif- 
ferent sizes, compared to the case of a static disk. However, as can 
be seen from Figure |5] although the eccentricities tend to group 
near eg for a wide range of s, the pericenters do not show such a 
general aligimient. Thus, it is difficult to say at this point whether 
the precession of the disk will help or hinder the accretional pro- 
cess. A more detailed analysis is thus necessary. 



5.2 Disk Density Profile 

It is important to note that in all our simulations presented so far, 
we have assumed a = 2 AU and a gas volumetric density of 
Pg = 5 X 10-^° gr/cm^ Since the gas drag scales linearly with the 
gas density pg, the drag dynamics of a given planetesimal scales 
with s/ Pg. For example, the largest limit cycle in Figure |5] corre- 
sponds to s = 0.5 km for the chosen gas density. However, if the 
gaseous disk becomes smaller by a factor A'^, the orbital motion 
of that planetesimal would become similar to one with a size of 
s = 0.5/N km. This suggests that in order to be able to assess the 
accretional probability of planetesimals, it is important to consider 
a realistic gas density profile for the gaseous disk. 

As shown by the results of hydrodynamical simulations of a 
gaseous disk around the primary star of the 7-Cephei system by 
Paardekooper et al. (2008) and Kley and Nelson (2008), one of the 
most important consequences of the secondary star is to truncate 
the original circumprimary disk at a point close to 5 AU from the 
primary and to introduce both an eccentricity and precession rate 
to the gaseous disk. Although the eccentricity and precession rates 
seem to depend on the disk parameters (Kley et al. 2008), the trun- 
cation radius and final density profile of the disk seem to be more 
robust and independent of the initial values. In particular, Kley and 
Nelson (2008) show that after ~ 100 orbital periods of the sec- 
ondary, the surface density profile acquires a near-equilibrium state 
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Table 1. Relative collision velocities (in m/s) for planetesimals with a = 2 
AU, Eg = 0.02 and for the precessing disk 2Tr/\gg\ = 1000 yrs. Radii ai'e 
in kilometers. 

which is practically linear with r, reaching a zero value at r ~ 5 
AU. 

In this paper, we adopt the conclusions by Kley and Nelson 
(2008), and for all subsequent applications of our simulations to 
the 7-Cephei system, we employ a gas density law of the form: 

E(r-) = B(ag(out) - Og) (33) 

where ag^out) = 5 AU is the outer edge of the disk and i? is a con- 
stant depending on the total mass of the disk. Since the gaseous 
disk is expected to be eccentric, its density should be constant 
along lines of constant semimajor axis ag and not along fixed val- 
ues of radial distance r. For this reason, we have expressed l l33t 
in terms of ag. However, it must be noted that the transformation 
E(r-) — >■ E(ag) makes use of our assumption that all the gas ele- 
ments share the same longitude of pericenter at any given time. In 
a more realistic gas disk this may not be the case. 

Also, because the density of the gaseous is not very steep near 
the origin, we can approximately express the total mass of the disk 
(Mt) as: 

Mt ~ jBag(out) (34) 

The value of the constant B can be explicitly calculated from this 
equation. The gas surface density E at 1 AU can now be written as 
Eo = -B(ag(out) ~ 1) which implies that the volume density pg is 
equal to 

In this equation, Hr = 0.05 is the scale height of the disk, 
taken constant for all values of ag. To specify the total mass of 
the gas disk, we note that the known planet in the binary system 
of 7-Cephei has a mass of ~ 1.6A/jup. We, therefore, choose 
Mt = 3Mjup. This value is larger than the one used by Kley 
and Nelson (2008) and smaller than the one by Paardekooper et al. 
(2008) and results in B = 1.3 x 10"" g/cm^. 

5.3 Encounter Velocities 

Tables[T]and[2]show values of encounter velocities between bodies 
of different sizes but with the same semimajor axis a = 2 AU. We 
compare two cases: a static disk (gg — 0) and a retrograde precess- 
ing disk with 27r/|(7g| = 1000 yrs. Table[T]corresponds to an almost 
circular gas eg = 0.02 while in Table|2]we considered a more ec- 
centric disk: eg = 0.2. The relative velocity at the encounter was 
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Table 2. Relative collision velocities (in m/s) for planetesimals with a = 2 
AU, Eg = 0.2 and for the precessing disk 2TT/\gg\ = 1000 yrs. Radii are 
in kilometers. 
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Figure 7. Relative collision velocities, as function of the semimajor axis, 
for several pairs of different sizes (values given in kilometers). Line colors 
correspond to different gas eccentricities, indicated in the upper left-hand 
plot. The broad horizontal brown line gives the limit for disruption colli- 
sions (Stewart and Leinhardt 2009). 

calculated using the equations of Whitmire et al. (1998). Since the 
precessing disk introduces a limit cycle in the orbits, the impact ve- 
locity AV will also oscillate. We can characterize this spread by its 
average (AV) and its minimum value AVmin- 

As shown in Table[T] encounter velocities are larger in a quasi- 
circular precessing disk, related to a larger eccentricity dispersion 
noted for s < 2 km (see Figure [5j. Such high encounter velocities 
hinder the accretion of small planetesimals in disks where ijg 7^ 0. 



In a static and more eccentric disk (Table|2j, the range of pos- 
sible equilibrium eccentricities lies between e/ and Cg (Figure|6]l. 
This means that, the larger the gas eccentricity, the broader the 
eccentricity dispersion. The encounter velocities in this case are 
higher and the collision among planetesimals may results in erosion 
and fragmentation. In an eccentric precessing disk, on the other 
hand, the disk precession reduces planetesimals' eccentricities to 
values below e / and causes bodies with different sizes to acquire 
similar eccentricities (bottom panels of Figure|6]l. In this case, even 
though Eg may be much higher than e/, the disk precession reduces 
the relative velocities during encounters even for small size bodies. 

As shown in Table |2] even with the eccentricity-damping ef- 
fect of the disk precession, the impact velocities for objects with 
radii s < 5 km in an eccentric disk are still larger than the dis- 
ruption velocity. However, it is important to note that these results 
were obtained for planetesimals at a = 2 AU. Because a truncated 
disk will have a much smaller gas density in its outer regions, it is 
possible that different semimajor axes will give different results. 

Figure|7]shows the variation of AV for six different planetes- 
imal pairs, as a function of the semimajor axis and for three values 
of Eg. In horizontal dashed lines we have also plotted the critical 
relative velocity for a catastrophic disruption V^o using the recipe 
developed by Stewart and Leinhardt (2009) for weak aggregates. 
As shown in this figure, AV decreases sharply in the outer parts of 
the disk for small planetesimals. Although for collisions between 
very small bodies the impact velocity is still above V^r, (upper left 
panel), we find that AV < V^jj for planetesimals with s ~ 2 
km and semimajor axis beyond 3 AU (upper right panel). Colli- 
sions between larger bodies are even more favorable (middle and 
lower panels). In all these cases the outcome of a collision seems 
to be accretion, at least in a significant portion of the disk. Even if 
the inner disk still appears hostile, it is possible to envision a sce- 
nario in which accretional collisions between small planetesimals 
occurs preferentially in the outer disk. Then, as these bodies grow 
and reach the inner regions due to orbital decay with the gas, they 
could continue their growth closer to the star. 

5.4 Timescales for Orbital Decay and Secular Equilibria 

In the previous section we analyzed the relative velocity of two 
planetesimals after they acquired their limit cycles. However, de- 
pending on the gas density and the radius of the object, the time 
required to reach the equilibrium solution may be longer than the 
typical collisional timescale (Paardekooper and Leinhardt 2010). 
During this time, the planetesimals may also undergo orbital de- 
cay. 

Although the interplay between dynamical and collisional 
evolution can only be evaluated with full numerical simulations, 
here we present an estimate of the characteristic timescales of the 
secular equilibria (re) and semimajor axis decay (r^). Starting with 
initial circular orbits and a = 2 AU, we define Te as the time nec- 
essary for a planetesimal with radius s to reach the final limit cycle 
with a relative error less than 1%. Since there are no equilibria in 
the semimajor axis, we define Ta as the time necessary for the body 
to decrease its semimajor axis by Aa — 1 AU. Figure [8] shows Te 
and Ta as functions of the planetesimal radius s, for different val- 
ues of the disk precession rate. The case of a non-precessing disk 
is identified with black lines, while color curves represent different 
precession rates (see figure caption for details). In the two top pan- 
els, we have assumed an almost circular gas disk (e^ = 0.02). As 
shown here, even a slow precession rate causes a significant reduc- 
tion in Te, which, for instance for a s = 0.1 km planetesimal, falls 
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Figure 8. Characteristic timescales for secular evolution in {k, h) variables 
(te) and orbital decay (Ta), as function of the planetesimal radius s, for 
several values of the disk precession rate. Top plots were drawn assuming 
a quasi-circular gas disk (eg = 0.02) while in the two bottom graphs we 
adopted eg = 0.2. In all cases the initial semimajor axes of the planetes- 
imals were chosen as a = 2 AU, and gas density given by equation )35t . 
Line colors are indicative of disk (retrograde) precession periods: 200 years 
(blue), 1000 years (green) and 4000 years (red). The case of a static (non- 
precessing) disk is shown in black. 
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Figure 9. Two orbital simulations of a planetesimal with s = 10 km 
in initially circular orbit with a = 2 AU. Black curves correspond to a 
static disk while red shows results using a retrograde precessing disk with 
'^'""/ISgl = 1000 yrs. In both cases the gas eccentricity was eg = 0.2. The 
continuous lines in the lower right-hand plot shows the center of the limit 
cycles as function of the semimajor, while the dashed lines show the max- 
imum and minimum values of the eccentricity. Notice the excellent agree- 
ment with the numerical simulation shown in the left-hand plot. 



from ~ 10^ years to ~ 10'^ years. Hence, a precession in the disk 
causes a much faster evolution from the initial conditions towards 
the limit cycle. 

Figure [8] also shows a comparison between Te and the 
timescale for orbital decay. For a static disk Te > r^, implying that 
the decay in semimajor axis occurs faster than the time necessary 
for the planetesimal to reach the secular solution. The opposite, 
however, occurs for a precessing disk, where now Te. < Ta for all 
values of s. In this case then, it is expected that the body will reach 
the limit cycle before suffering any significant orbital decay. 

Figure[9]shows the results of two sets of N-body simulations. 
Planetesimals were initially at a = 2 AU with e = and an 
adopted radius of s = 10 km. The gas disk was assumed to have an 
eccentricity of = 0.2. Black curves correspond to a static disk 
(no precession) while red curves denote a disk with (retrograde) 
precession period of 1000 years. As expected from Figure [S] the 
semimajor axis of a planetesimal falls more rapidly for a precess- 
ing disk, reaching 1 AU in timescales slightly over 10^ years. The 
orbital decay rate increases with time due to the greater gas den- 
sity near the central star (see equation (135b). The evolution of the 
eccentricity in this case shows different behaviors. While the plan- 
etesimal approaches a quasi-circular orbit in the precessing disk, 
the opposite occurs in the static case, reaching values close to Cg 
for small values of a. 

The two bottom panels of Figure |9] show the eccentricity as 
a function of the semimajor axis. On the left-hand panels we have 
again plotted the results of the exact N-body simulations. Due to the 
orbital decay, the orbital evolution occurs from the right to the left 
of the graph. On the right, we plot the values of the center of the 
limit cycles (continuous lines) as well as the minimum and max- 
imum values of the eccentricity (dashed lines), each determined 



numerically from the averaged model ([28} for fixed values of the 
semimajor axis. 

Except for the first few tenths of AU close to the initial con- 
dition, the rest of the planetesimal's evolution occurs very close to 
the instantaneous limit cycles for each value of the semimajor axis. 
In other words, even in the presence of the significant orbital decay, 
the secular dynamics of the planetesimals is expected to be dictated 
by the limit cycles in (fc, h). Moreover, except for the first few 10^ 
years, it is expected that the encounter velocities of collisions be- 
tween planetesimals be determined by the equilibrium solutions of 
the secular problem, and not by the transient evolution from the 
initial conditions to the limit cycles. 



6 CONCLUSIONS 

In this paper we have presented an initial study of the dynamics 
of individual planetesimals in circumstellar orbits of a tight binary 
system. Apart from the gravitational perturbations of the secondary 
star, we have also considered the effects of a non-linear gas drag 
due to an eccentric precessing gaseous disk. 

A non-zero precession frequency gg in the gas disk introduces 
three main changes in the secular dynamics of the solid bodies. 
On one hand, the eccentricity and longitude of pericenter no longer 
reach stable fixed points (in the averaged system), but display peri- 
odic orbits in the (fc, h) plane with the same frequency as the disk. 
This appears to be independent of the semimajor axis of the plan- 
etesimal, and thus of the secular frequency g of the gravitational 
perturbations due to the secondary star. The amplitude of the limit 
cycle is inversely proportional to the planetesimal radius s, reach- 
ing values close to zero for large bodies. 

Another issue is the range of planetesimals' final eccentrici- 
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ties. For large planetesimals, the limit e — > e / is still observed for 
all values of gg, where e/ is the forced eccentricity induced by the 
secondary star. However, the behavior of small bodies is dependent 
on Qg.ln a static disk, e — > eg for planetesimals strongly coupled 
to the gas. However, the final eccentricity is significantly lower in 
a precessing disk, reaching almost circular orbits for high values of 
9g- 

Finally, we have also noted that the precession also reduces the 
time necessary for a given planetesimal to reach the stationary limit 
cycle, starting from initially circular orbits. Although the orbital de- 
cay is also accelerated, the characteristics timescale associated with 
the semimajor axis is always longer than the secular timescale. This 
seems to imply that, except for a very small time interval immedi- 
ately following the initial conditions, the full dynamics of the plan- 
etesimals should be well represented by the equilibrium periodic 
orbits obtained for each instantaneous value of a. 

Applying these results to the 7-Cephei binary system, and as- 
suming a gas density profile similar to those obtained from hydro- 
dynamical simulations, we find that the precession of the disk also 
causes significant differences in the relative velocity of colliding 
planetesimals of different radii. For quasi-circular gas disks, a non- 
zero value of Qg appears to increase the encounter velocities, hin- 
dering the possibility of accretional collisions among small ~ 1 — 5 
km bodies. However, the opposite effect is observed for values of 
Cg ~ 0.2, where a non-zero precession of the gas actually helps 
reduce the eccentricity dispersion and reduce the coUisional veloc- 
ities. 

Assuming a gas disk of total mass ~ 3Af,jup we found that 
accretional collisions may occur in the outer parts of the disk for 
s > 2 km. For larger bodies, the location of minimum relative ve- 
locities appears closer to the central star. However, full N-body sim- 
ulations are necessary in order to fully evaluate the consequences 
of this dynamics on a mutually interacting protoplanetary swarm. 
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